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Significant microstructural anisotropy is known to develop during shearing flow of attractive particle suspensions. These sus¬ 
pensions, and their capacity to form conductive networks, play a key role in flow-battery technology, among other applications. 
Herein, we present and test an analytical model for the tensorial conductivity of attractive particle suspensions. The model uti¬ 
lizes the mean fabric of the network to characterize the structure, and the relationship to the conductivity is inspired by a lattice 
argument. We test the accuracy of our model against a large number of computer-generated suspension networks, based on 
multiple in-house generation protocols, giving rise to particle networks that emulate the physical system. The model is shown to 
adequately capture the tensorial conductivity, both in terms of its invariants and its mean directionality. 


Introduction 

The electrical conductivity of heterogeneous materials has 
been extensively studied by many different researchers over 
the years 34 12 1920 . The literature primarily focuses on het¬ 
erogeneous materials which are mixtures of two materials 
that each have different, isotropic electrical conductivities. 
The most well-known result is that of Maxwell, which is 
based on an effective-medium approximation for dilute sus¬ 
pensions 14 . Hashin and Shtrikman approached the problem 
in a different way. Rather than attempt to solve for an exact 
expression for the effective conductivity of a randomly struc¬ 
tured material, they applied a variational method to derive up¬ 
per and lower bounds on the effective conductivity 9 . They 
chose to use a variational approach to derive bounds on the 
conductivity because solving the exact problem for an arbi¬ 
trarily structured heterogeneous material was analytically in¬ 
tractable. Torquato 19 21 has studied the effective conductivity 
problem in great depth. He has improved the bounds laid out 
by Hashin and Shtrikman, has solved for effective conductiv¬ 
ity of a number of different lattice types, and has expressed 
the exact tensorial effective conductivity in terms of an infinite 
series of N-point probability functions, which can be used to 
describe the microstructure of a heterogeneous material. The 
particular case of a suspension consisting of a conductive par¬ 
ticle network within an insulating medium has been consid¬ 
ered theoretically, to our knowledge, in one existing study 12 . 
The approach they take assumes a spatially homogeneous po¬ 
tential gradient field imposed upon the structure, leading to a 
model for the conductivity that can be proven to be an upper 
bound. 

Much of the aforementioned work is concerned with the 
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isotropic conductivity of heterogeneous materials. In this 
work, we aim to model the full tensorial conductivity, with a 
focus on suspended networks of conductive particles. These 
particle networks are of practical importance, especially in 
flowable battery technology currently under development by 
the Joint Center for Energy Storage Research (JCESR) 5 . In 
these batteries, a conductive, flowing suspension of carbon 
black forms an integral component of the system, see Fig¬ 
ure [lj a). It has been shown in related systems 11 that shear¬ 
ing flows induce anisotropy in a contact network of sus¬ 
pended particles, as pictured in Figure |TJb). In instances 
where suspension conductivity arises from particle-particle 
contacts, this structure anisotropy should give rise to conduc¬ 
tivity anisotropy. It is this behavior that we seek to describe. It 
has been shown experimentally that the electrical conductivity 
of a suspension is highly sensitive to shear rate 1 , dropping by 
several orders of magnitude as shear rate increases. From this 
observation and the evidence of particle microstructure chang¬ 
ing in shearing flow, we deduce that a suitably chosen descrip¬ 
tion of the particle network should be sufficient to predict the 
electrical conductivity of a suspension. 

In the granular media literature, a great deal of attention has 
been given to describing the structure of the contact network 
between particles. Perhaps the simplest structural measure 
for such a network that includes anisotropy is the fabric ten- 
sor 151618 . While more complex structural measures exist, 
such as pair- and higher-order particle correlation functions^- 19 ! 
whose use could enable greater accuracy in constructing a 
conductivity model, we shall show that a suitable model can 
be achieved solely in terms of the fabric. Key to our model de¬ 
velopment is the solution of a simple case, based on a network 
conforming to a lattice structure. The results instruct the form 
for a new conductivity model, whose accuracy is then tested 
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Fig. 1 (a) Image of a carbon black particle network, an elec¬ 
trically conductive suspension®!, (b) Image of an effective two- 
dimenionsional suspension (attractive polystyrene beads on a fluid 
surface), which has been subjected to shearing. Note the formation 
of an anisotropic contact network between particles.EI 


against many thousands of random particle networks. To ex¬ 
plore a range of particle networks, we describe two distinct 
algorithms for creating random packings — one for denser 
packings, and one for more dilute packings that closely re¬ 
semble those formed by carbon-black — and demonstrate the 
model’s predictive capability against thousands of packings 
generated from both algorithms. 


Homogenization 


The tensorial form of Ohm’s law relates the electric field vec¬ 
tor E to the current density vector J through a second-order 
conductivity tensor K, i.e. 


J = KE (1) 

The conductivity tensor is a symmetric, positive-definite ten¬ 
sor 21 . An effective conductivity for a representative volume 
O of a heterogeneous material must be defined prior to any 
analytical or numerical work. The effective conductivity of an 


ergodic medium is defined by 


(J) = K(E) 

where (E) and (J) are, respectively, the spatially-averaged 
electric and current density fields over To avoid a 

possibly over-reaching assumption of ergodicity — our tests 
will be conducted on finite domains — we specify that 
(E) is imposed by prescribing a linear boundary potential 
(p(x £ dCi) = —(E) • x, and that (J) is redefined as the flux that 
is power-conjugate to (E). That is, 

(E)-(J) = i|-V<p.jdV (3) 

a 

where j is the local current density field. In the ergodic limit of 
the ensuing analysis, (J) reduces to a standard spatial average. 

Assuming that the current density obeys Kirchoff’s current 
law and Ohm’s law — respectively, V • j = 0 and j = —aV(p 
for some non-negative conductivity field cr(x) —a symmetric, 
positive-definite conductivity tensor K must exist that obeys 
([2]). By using a calculus identity, Eq[3]can be transformed into 

(E)-K(E) = i|-<pj.ndA (4) 

dQ. 


where n is the outward-pointing normal vector. 


We model the particles as perfect conductors, the fluid as a 
perfect insulator, and we suppose electrical resistance arises 
only at the contacts between particles. Likewise, the field <p is 
approximated as a constant within each particle but possibly 
varying from particle to particle. The above integral can now 
be broken into a sum of integrals over the boundary. In the lo¬ 
cations where the boundary passes through free space (i.e., not 
a particle), then we know that j is exactly 0. This leaves only 
the parts of the boundary that pass through particles, which al¬ 
lows us to write the integral over the set of boundary particles 
B, i.e. 


(E)-K(E) 


1 

V 


L-w- 


/ 


j 


ndA. 


(5) 


where fi, is the intersection of the ;th boundary particle with 
d £2, and the potential within particle i, denoted <p, above, can 
be brought outside the integral since it is constant within a 
particle. Although the precise nature of j is unknown within 
the particle, the value of the integral f dn j - ndA is the current 
that is flowing out of £2. Denoting this current as I° m we can 
write the final expression for the right-hand-side of 


<E>-K(E) = i£-<p,/r (6) 

v ieB 
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Fig. 2 Schematic of a small resistor network with nodes and edges 
labeled according to our conventions. 



Fig. 3 Schematic of particles in contact showing contact vectors n,. 


The three independent components of K,, can be determined 
by performing multiple simulations on the same particle net¬ 
work with three non-colinear choices of (E). 

By our assumptions for the particle properties, the problem 
can be reduced further to that of a resistor network. The net¬ 
work is defined by the set of particles acting as the nodes, 
which are connected by a set of contacts acting as the edges, 
which carry a resistance R c . A schematic of an example net¬ 
work with 9 nodes and 12 edges can be found in figure [2] 
Supposing an /V-particle sample and letting i, n n represent the 
(signed) current flowing from particle m to n, Ohm’s and Kir- 
choff’s law can be rewritten in their simpler discrete form. 


_ *Pm tyn 

Rc 


and 


y\im,n = 0 for all m. 

n 


(7) 

( 8 ) 


can be defined for particle p by the relation'! 5 — 

A p = }Tn/®n,- (9) 

i 

where 0 denotes the dyadic product, and n, is the unit normal 
vector connecting particle centroids of the i’th contact on the 
particle. This is illustrated in Figure [3] To homogenize over 
the entire particle network, or at least meso-sized region of it, 
the average fabric tensor is defined as the system average of 
the particle fabric tensors. 


A = 


1 


Aparticles 


^particles 

E A p 


p= 1 


( 10 ) 


The definition of the fabric tensor has some attractive features. 
It is symmetric and positive-semidefinite, guaranteeing that 
the eigenvalues are non-negative and that the eigenvectors are 
orthogonal. These properties are shared by the conductivity 
tensor K, suggesting the fabric tensor could be an appropriate 
independent variable in the conductivity’s functional form. 


Lattice-Reduced Model 


We propose an analytical model to elucidate the connection 
between electrical conductivity and the fabric tensor based on 
a simplified lattice structure. We will test this model’s appli¬ 
cability to random packings in the later sections. 

The particles are imagined to live on an idealized infinite, pe¬ 
riodic lattice. The lattice is parameterized by a set of numbers 
that describe the particle size and spacing. These parameters 
are (1) particle diameter D p , (2) distance in x-direction be¬ 
tween chains d x , (3) distance in y-direction between chains 
d y , (4) distance in z-direction between chains d z . In 2D, only 
the first three parameters are used. An illustration of a 2D lat¬ 
tice characterized by these parameters is shown in figure[4|a), 
with its fundamental unit cell shown in figure[4jb). 

Both the average fabric tensor and effective conductivity can 
be computed analytically. The average fabric tensor is defined 
as the spatial average of the fabric tensor for all of the particles 
in the unit cell and ultimately results in the formua 


A = 


2 

N.X +Ny — 1 


N x 0 
0 Ny 


(ID 


Solving these linear equations for a given particle network en¬ 
ables us to calculate I out in (|6]l and hence the conductivity ten¬ 
sor for the network. 

We choose to use the fabric tensor as the measure of the net¬ 
work structure. The particle-level fabric is a local quantity that 


In this expression, the key quantities to recognize are the num¬ 
ber of particles in the x-oriented chain, N x = d x /D p , and the 
number of particles in the y-oriented chain, N y = d y /D p . 

Next, the effective conductivity was derived for the unit cell. 
To do this, imagine applying an arbitrary voltage difference 
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Fig. 4 Idealized particle lattice and unit cell from which fabric- 
conductivity relation was derived, (a) Example 2D idealized parti¬ 
cle lattice, (b) 2D lattice unit cell and its resistor network analog. 
Neighboring unit cells are shown in gray dashed lines. 


across the x-oriented and y-oriented chains separately. These 
voltages are A (p x and A<p v , respectively, effective resistance. 
By applying Ohm’s law through the corresponding chains, we 
can recover the components of the vector form of Ohm’s Law 
shown in (JTJ). For example, for the x-oriented chain 

jx= Gp;) ( v<p)a ' (12) 

with (V<p) x = A(p x /d x . Due to the geometry of the problem, 
we know that the off-diagonal components of the conductivity 
tensor K are exactly zero. Therefore, we can say 


*n = 


1 

NyRc' 


Similarly analysis yields 


*22 


1 

NxRt 


(13) 


(14) 


Finally, the parameters N x and N y can be algebraically elimi¬ 
nated to give the components of K in terms of the components 
of A, yielding the tensorial relationship 


K = 


1 trA - 2 

_ A 

R c detA 


(15) 


We refer to the formula in GD as the “lattice model”. A simi¬ 
lar analysis can be carried out for a three-dimensional unit cell. 


which will yield the following expression for the conductivity 
tensor. 


K = 


1 (trA —2) : 
4 D p R c detA 


A. 


(16) 


The formulae above apply when trA — 2 is non-negative. Oth¬ 
erwise the solution is K = 0. 


Despite its inspiration from the lattice structure, there are sev¬ 
eral reasons to consider the applicability of the lattice model to 
more general particle networks. For one, the formula purports 
codirectionality of the fabric and conductivity, i.e. the devia¬ 
tors of the two tensors are aligned, implying that the direction 
of anisotropy of one tensor gives the anisotropy direction of 
the other, which to a first approximation ought to match the 
behavior of general particle networks. Second, the results im¬ 
ply that conductivity should vanish when trA < 2, which is 
sensible more generally (though not strictly) because particles 
in a percolating chain, as needed to conduct current across the 
sample, must have coordination number at least two. Above 
this threshold, conductivity increases with trA in line with 
one’s basic intuition for more highly coordinated networks. 

We are aware of one other fabric-based analytical model for 
conductive particle networks, which was developed by Jagota 
and HuiEl. In their work, a uniformity hypothesis is made with 
regard to the potential gradient, which results in a conductivity 
model that is fully linear in the fabric tensor. 


K = 



(17) 


The above, which can be proven to be an upper-bound on the 
real conductivity, is for a two-dimensional system and Ny is 
the particle number fraction (per area in 2D). This model dif¬ 
fers from ours most notably in that the conductivity is not 
thresholded by the coordination number, the formula depends 
explicitly on the particle area fraction as well fabric, and it 
does not depend on the fabric determinant. 


Numerical Simulation 

In order to perform numerical experiments and determine the 
generality of the lattice model, a large number of random par¬ 
ticle networks (packings) must be created. There are a num¬ 
ber of methods to do this already in the granular and partic¬ 
ulate matter literature. See the references for a broad sum¬ 
mary of the currently available granular packing algorithms 2 . 
Attractive suspensions have been modeled with the Diffusion- 
Limited Aggregation (DLA) model of Witten and Sander 13 . 
A common feature of many of the granular statics methods is 
that they solve force equilibrium equations for a system of par¬ 
ticles This was not a feature that was required for this study. 
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Fig. 5 Example (using a small number of particles) of a dense particle 
packing resulting from Algorithm^!] 

so these types of methods were not used, in the interest of 
saving computational time. Instead, we developed two meth¬ 
ods for creating two-dimensional random contact networks of 
particles, and we tested our model against numerous packings 
generated by each method. Both methods allow us to influence 
the resulting anisotropic structure of the packings. 

Algorithm 1: Our first packing algorithm was designed to cre¬ 
ate a dense random contact networks of particles. This is in 
contrast to a later algorithm, to be described below, which cre¬ 
ated packings that resulted in much lower-density packings. 
The dense packings were created by perturbing a 2D hexago¬ 
nal close-packing of particles. This was achieved by placing 
points into a triangular lattice, adding random noise to the po¬ 
sition of each point, and finally growing each particle as large 
as possible such that no particles overlapped. Anisotropy can 
be influenced by shearing the points with an affine transfor¬ 
mation x' =Fx before growing the radii. This process is de¬ 
scribed in pseudocode below (Algorithm [TJ. An example of 
the resulting packing overlaid by its analogous resistor net¬ 
work is shown in figure [5] 


Algorithm 1 

Seed Lx L box with close-packed points 
Perturb points with random noise 
Move each point to new location x' by x' = Fx 
while Not all radii frozen do 

Find smallest distance that any particle can grow 
Grow all particles by this amount 
Freeze radii of particles that come into contact 
end while 


Algorithm 2: This procedure was motivated by a need to better 
understand the conductivity of carbon black suspensions in an 


insulating medium. The self-attraction carbon black particles 
leads to fractal particle networks that are electrically percolat¬ 
ing at low volume fraction (below 1 vol%)®. 

To produce structures that more closely resemble carbon black 
suspensions, we developed our second packing algorithm, 
which is inspired by the “hit-and-stick” behavior of the car¬ 
bon particles. In addition, the new algorithm is able to include 
the effects of particle Brownian motion but this is not essential 
to the algorithm. 

First, clusters (single particles at this stage) are seeded ran¬ 
domly into a L d box, where d is the number of spatial dimen¬ 
sions. Next, a linear velocity field is imposed directly on each 
cluster’s centroid according to 

v = -B(x-0) (18) 

where O is a point in the middle of the original box. This im¬ 
posed velocity field serves to pull all of the clusters together. 
The matrix B is a d x d matrix that allows us to impose an 
anisotropic velocity field. This allows us to influence (but 
not completely impose) the fabric tensor that results from this 
packing method. After the velocity field is imposed, the par¬ 
ticle positions are updated by assuming a time step dt (com¬ 
puted at runtime). Then, the clusters are checked to determine 
whether any contacts have been made with other clusters. If 
so, the clusters are cohered into a single cluster for all future 
steps. This process of imposing velocity, updating positions, 
and handling contacts is repeated until only a single cluster re¬ 
mains. The process is outlined in pseudocode in Algorithm[2] 
An example of a packing resulting from this process is shown 
in figure[6]and a larger example is displayed in figure[7] 

The box-counting fractal dimensioned'the resulting packings 
was computed in order to determine if they resembled real- 
life packings found in experiments. The fractal dimension of 
packings produced by this method is approximately d = 1.75. 
This was compared against the particle network image in fig¬ 
ure |T] This network has a fractal dimension of approximately 
d = 1.7 ±0.1. Uncertainty in the measurement is due to the 
image processing techniques used to identify particles. Based 
on these measurements, we are satisfied that this algorithm 
produces realistic packings, although more detailed correla¬ 
tion function measurements would be needed for a firmer con¬ 
clusion. 

Applying boundary conditions: In order to apply the solution 
method described above to an arbitrary packing of particles, 
appropriate boundary conditions must be applied. In these 
simulations, a prescribed voltage was applied to particles all 
around the boundary. This process consists of two steps: first, 
the boundary must be identified, and second, the linear system 
must be updated to reflect the known voltages. 
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Algorithm 2 

Seed N clusters (particles) in L d square 
While NClusters > 1 do 

Move clusters according to v = — B(x — O) 

Locate collisions between clusters 
Combine clusters in contact and recompute centroids 
end while 



Fig. 6 Example (using a small number of particles) of a packing 
resulting from Algorithm[2]using a small number of particles. 


For the first packing algorithm, identifying the boundary is a 
trivial process, since the particle locations are known a priori . 
For algorithm[2] however, the particle positions are not known. 
A boundary can be located visually quite easily at the end 
of the simulation process, but performing this step manually 
would be prohibitively slow. In order to expedite and automate 
the simulation process, the following method was devised to 
locate the boundary. 

First, histograms of the particle x and y positions were sepa¬ 
rately created. To find the “left” and “right” boundaries, de¬ 
noted x~ and x + respectively, the histogram of x positions was 
thresholded. The value x is defined as the smallest x value 
where the histogram reaches 50% of its maximum value. The 
value x + is defined as the largest x value that meets the same 
criterion. The top and bottom boundaries, y + and y~, are 
found in the same manner using the histogram of particle y 
coordinates. The threshold value 50% was determined emper- 
ically to locate the same boundary that one would identify vi¬ 
sually. An example packing and its associated x-position his¬ 
togram is shown below in figure [7] to demonstrate the efficacy 
of the method. Once the location of the boundary has been 
identified, all particles whose centers fall less than one radius 
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Fig. 7 Example 10,000-particle packing (from Algorithm 2) with 
its associated x-position histogram and the boundary selected by the 
method. 


away from the lines are marked as being “boundary particles”. 

The expression in ([6]) can be computed easily from the solu¬ 
tion of the particle network, so by judiciously choosing (E), 
the components of K,, can be extracted. In two dimensions, 
the effective conductivity tensor has three independent com¬ 
ponents, so three simulations are sufficient to extract all of the 
components. The K\\ component can be extracted by setting 
(E) = e T . This corresponds to evaluating the integral for an 
applied boundary voltage of (p = —x. The remaining tensor 
components may be similarly extracted by applying specific 
potential fields at the boundary and evaluating the summation 
given in •Ill- 


Tests 


The previously described packing algorithms and solution 
procedures for the current/potential have been implemented 
in Matlab. Algorithm [I] was used to create 50,000 separate 
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400-particle packings. In all of these packings, the F\ \ and 
F 22 components of the affine transformation F equalled 1.0. 
The F \2 component that controlled the shearing of the packing 
ranged between 0 and 0.5 in increments of 0.01. Any par¬ 
ticles that were sheared out of the original bounding rectan¬ 
gle were reflected to the other side of the box to return the 
packing to a rectangular geometry. Algorithm [2] was used to 
create 10,000 separate 5,000-particle packings. In the B ma¬ 
trix, the B 11 component remained 1.0, and the B 22 component 
was varied in [1.0,1.9] in increments of 0.1 to influence the 
level of anisotropy of the resulting packings. In all simula¬ 
tions, the contact resistance R c was assigned to be 1, so it did 
not have any affect on the following analyses. After applying 
the previously described procedure to each packing to obtain 
the effective conductivity tensor and average fabric tensor for 
each packing, the data were analyzed to determine how well 
the results agree with the model’s prediction for the isotropic 
magnitude, the deviatoric magnitude, and the direction of con¬ 
ductivity. These tests are described next, and thereafter we 
shall proceed to show how well the lattice model performs 
compared to the existing model, equation ( fl7] >. 

The isotropic behavior of the conductivity can be investigated 
by taking the trace of both sides of ( fT5| . The average coordi¬ 
nation number is the most natural independent variable when 
examining the isotropic behavior, so in addition to taking the 
trace of both sides of ( fl5| ), both sides were multiplied by det A 
in order to make the right-hand side a single-valued function 
of trA. This results in (fl9|. 


R c trK det A = (trA — 2) trA (19) 

The results of the simulations are plotted together with the an¬ 
alytical curve given by ( fl9| in figure [8] It was found that the 
analytical solution is usually an upper bound on the measured 
conductivity. This can be explained by the fact that the ana¬ 
lytical model was derived from an idealized system where the 
chains span a unit cell in a straight line. Since the total resis¬ 
tance of a chain is proportional to the number of contacts in the 
chain, it follows that the shortest chain between any two points 
is the lowest resistance path, and therefore most conductive. 
Since the model was derived from a straight-chain idealiza¬ 
tion, it implies an upper bound on the conductivity. This logic 
is less valid in low-coordinated systems, which have many dis¬ 
connected groupings of one or two particles; low-coordinated 
systems rarely if ever occur from Algorithm 2 or in actual car¬ 
bon black suspension networks. In this case, the trace of the 
system’s fabric can be less than 2 but percolating chains may 
still exist to produce small but non-zero conductivity. This 
effect is evident in the figure in the data of Algorithm 1. 

Next, we determine the extent the analytical lattice model pre¬ 
dicts the anisotropy of the conductivity. To remove the influ¬ 



Fig. 8 Predicted relationship between the (modified) trace of the 
conductivity and the fabric trace, compared to numerical results of 
50,000 packings generated by Algorithm 1 and 10,000 generated un¬ 
der Algorithm 2. Inset is a zoom-in of the vicinity of trA = 2. 


ence of the isotropic behavior, we take the deviator of both 
sides of ( fl5} . In this case, the most natural independent vari¬ 
able is the magnitude of the fabric deviator, so the resulting 
equation was manipulated to be a single-valued function of 
this quantity. After manipulation, © can be written as ([20}. 


Rc Kq 


Ao / detA \ 
|A 0 | \trA — 2/ 


|Ao| 


( 20 ) 


where a subscript 0 denotes the deviator of the tensor, and the 
term is commonly referred to as the direction or sign of 
the tensor Ao- The left hand side of this was plotted against 
|Ao| to test the predictive power of the model. It can be seen 
in figure [9] that, although there is a large amount of noise in 
the measurements, the model captures the mean behavior very 
closely. 

The final prediction that must be examined is the notion of 
codirectionality. The analytical model in ( [15} predicts that 
the fabric and conductivity tensors have the same eigenvec¬ 
tors. To examine this, the angle difference between the fabric 
and conductivity deviators was calculated, which is equivalent 
to the (signed) angle between the eigenvectors corresponding 
to the largest eigenvalues of the two tensors, denoted and 
e A . The deviators were chosen because, in 2D, the eigenvec¬ 
tor corresponding to the positive eigenvalue can be unambigu¬ 
ously chosen. The probability density function of the angle 
difference as a function of A0 is plotted in figure [TO] It can 
be seen that this distribution is symmetrically centered around 
zero, indicating that the fabric and conductivity are strongly 
codirectional. 
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Fig. 9 Predicted relationship between effective magnitude of the 
anisotropy of the conductivity and the invariants of the fabric. Er¬ 
ror bars show ± one standard deviation. 



AO (radians) 

Fig. 10 PDF of angle differences are distributed around zero, in¬ 
dicating codirectionality of the fabric and conductivity tensors, as 
predicted by the analytical model, i.e. G3- 


Finally, we also compared the lattice model, ( fl5] l, to the ex¬ 
isting model by Jagota & Hui 12j shown in ( fT7j ). For the same 
60,000 packings generated using both packing algorithms, we 
computed the relative error of the prediction of the trace and 
the determinant of the conductivity using each model and 
plotted the results in figures 0 and [12] In every case, we 
found that the new lattice model predictions were closer to 
the true values from the numerical experiments than the pre¬ 
vious model by Jagota & Hui. On the other hand, the Jagota 
& Hui model maintains a strong upper bound on both invari¬ 
ants of the conductivity tensor, whereas the lattice model is 
not strictly an upper bound, as previously discussed. 




Fig. 11 Plot of the relative error of the trace of conductivity. (Left) 
Relative error from packings created with Algorithm [T] (Right) Rel¬ 
ative error from packings created with Algorithm]!] 



Fig. 12 Plot of the relative error of the determinant of conductivity. 
(Left) Relative error from packings created with Algorithm]!] (Right) 
Relative error from packings created with Algorithm]!] 


Discussion and Conclusions 


In this paper we have derived and tested a new model relating 
the structure of a packing of particles to its tensorial electrical 
conductivity. The assumptions implicit in the model are that 
the suspending medium is a perfect insulator and that electri¬ 
cal resistance arises only at particle contacts. The structural 
measurement used was the fabric tensor, and the model arises 
from a straightforward analysis of a representative problem 
involving a lattice structure. The resulting model takes a non¬ 
linear functional form, and was tested multiple ways against 
numerical simulations of many thousands of random parti¬ 
cle packings. The agreement in its predictions of the various 
scalar properties and tensorial orientation is significant, espe¬ 
cially in light of the simplistic nature of the fabric tensor being 
the sole independent variable for the model. In our tests, the 
lattice model’s accuracy was shown to be higher than an ex¬ 
isting conductivity model, a model which requires more struc¬ 
tural input data than the lattice model. While it is definitely 
possible to write a more accurate model by including depen¬ 
dences on more structural variables — some of our data spread 
is due to the finite nature of the datasets, but some is surely due 
to modeling error — the current simplicity of the lattice model 
is an advantage for its usage in engineering applications in¬ 
volving flowing suspension networks. Modeling frameworks 
for the evolution of anisotropy tensors in flowing media have 
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been developed over the last decades 7 M3 keeping our model 
in terms of fabric, then, suggests a path to the simulation of 
simultaneous flow and current transfer fields in nontrivial sys¬ 
tems by coupling a fabric evolution rule and a rheology with 
our conductivity model. Such a capability would be key in the 
targeted application of modeling flow battery systems, which 
rely on a flowing conductive suspension that closely resembles 
the idealized system that we considered. 
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